## This program generates Figure 8 from "Financial Constraints, Sectoral Heterogeneity,
##  and the Cyclicality of Investment" by Cooper Howes

##### Step 0: Preliminaries #####

## Clear variables and data
rm(list=ls())

## Load packages
library(tidyverse)

##### Step 1: Load data #####

## Load impulse responses from model (created in Dynare file "baseline_model.mod")
irfdata <- t(read.csv("irfs/shock_irfs_baseline.csv",header=F)[1:9,])

## Extract names and relabel
colnames(irfdata) <- irfdata[1,]
rownames(irfdata) <- NULL
irfdata <- irfdata[-1,]

## Convert IRFs to numeric
irfdata <- apply(irfdata, 2, as.numeric)

## Pick which impulse responses to plot
## Names are "variable"_"shock", with the names of each defined in the Dynare file
varnames <- c("k_d_e_m","k_n_e_m","k_e_m",
              "i_d_e_m","i_n_e_m","c_e_m",
              "p_d_e_m","pi_e_m","y_e_m")

## Full variable names to be used in plots
longvarnames <- c("Capital stock (durable)","Capital stock (nondurable)","Capital Stock (total)",
                  "Investment (durable)"," Investment (nondurable)","Consumption",
                  "Relative durable price","Inflation","Total output")

##### Step 2: Plot IRFs #####

## Plot everything in one big figure
windows(width=30,height=22)
par(oma=c(3,3,3,0))
numplots <- length(varnames)
m <- matrix(1:length(varnames),nrow = 3,byrow=T)
layout(mat = m)

qtrs <- 0:16

for(i in varnames){
  
  ## Extract the IRF to plot 
  ## Note that the data are transposed when loading and are labeled by row
  irf <- irfdata[,i] 
  
  ## Calculate y-axis range
  y_range <- range(irf)
  
  ## Plot each IRF
  par(mar = c(1.75,3.5,2,2))
  plot(x=qtrs,y=irf,type="l",col="blue",lwd=3,ylim=y_range,ann=FALSE,cex.main=1.25,
       axes=F,frame=T)
  axis(1,cex.axis=2)
  axis(2,cex.axis=2,at=pretty(y_range,n=5),lab=paste0(pretty(y_range,n=5)*100, "%"),las=TRUE)
  abline(a=0,b=0,col="black",lwd=1)
  title(main=longvarnames[which(varnames==i)], col.main="blue", font.main=2,cex.main=2)
  
  ## Manually add x-axis label to bottom-middle figure (inflation)
  if(i=="pi_e_m"){mtext(text="Quarters after shock",cex=1.5,side=1,line=3)}
  
  
}

## Add legend
par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
title(main="Impulse response to 100bp contractionary monetary shock",cex.main=2,col.main="blue",line=-2)

